Approche variât ionnelle pour le calcul bayésien dans les 
problèmes inverses en imagerie 
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Résumé — Dans une approche bayésienne non supervisée pour la résolution d'un problème inverse, on cherche à estimer conjoin- 
tement la grandeur inconnue / et les paramètres G. Ceci se fait en utilisant la loi a posteriori conjointe p(f, 0\g). L'expression de 
cette loi est souvent complexe et son exploration et le calcul des estimateurs bayésiens nécessitent soit l'optimisation des critères 
souvent non convexes ou le calcul d'espérances des lois non gaussiennes multivariées. Dans tous ces cas, il y a souvent besoin 
■ de faire des approximations. Nous avions déjà exploré les possibilités de l'approximation de Laplace et les méthodes d'échan- 
tillonnage MCMC. Ici, nous explorons l'approximation de p(f, 0\g) par une loi séparable en / et en G. Ceci permet de proposer 
des algorithmes itératifs plus abordables en coût de calcul, surtout, si on choisit ces lois approchantes dans des familles des lois 
exponentielles. Le principal objet de ce papier est de fournir des détails des différents algorithmes que l'on obtient pour différents 
choix de ces familles. 

Abstract — In a non supervised Bayesian estimation approach for inverse problems in imaging Systems, one tries to estimate 
jointly the unknown image pixels / and the hyperparameters G. This is, in gênerai, done through the joint posterior law 
p(f,0\g). The expression of this joint law is often very complex and its exploration through sampling and computation of the 
point estimators such as MAP and posterior means need either optimization of non convex criteria or intégration of non Gaussian 
and multi variate probability laws. In any of thèse cases, we need to do approximations. We had explored before the possibilities 
of Laplace approximation and sampling by MCMC. In this paper, we explore the possibility of approximating this joint law by 
a séparable one in / and in G. This gives the possibility of developing itérative algorithms with more reasonable computational 
cost, in particular, if the approximating laws are choosed in the exponential conjugate families. The main objective of this paper 
is to give détails of différent algorithms we obtain with différent choices of thèse families. 
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1 Introduction 

Dans une approche estimation bayésienne non supervi- 
sée pour résoudre un problème inverse, on commence par 
écrire l'expression de la loi a posteriori conjointe des in- 
connues / et des hyper-paramètres 8 : 



> 



Aussi, notant par 



et par 



Kh(q : p) 



p(f,8\g;M) = 



p(g,f,0\M) _p(g\f,0;M) P (f,0\M) 



P (g\M) 



p{g\M) 



(1) 

Dans cette relation p(g\f, 8; A4) est la vraisemblance des 
inconnues dont l'expression s'obtient à partir d'un modèle 
liant les inconnues aux données g (modélisation du pro- 
blème directe), p(f, 8\A4) est la loi a priori des inconnues 
et 

p(g\M) = H p(g\f, 8; M) P (f\0; M) p{8\M) df d8 

(2) 

est ce qu'on appelle l'évidence du modèle A4. 

Il est intéressant de mentionner que, pour n'importe 
quelle loi de probabilité q(f, 8) on a 



on montre facilement que 

p(g\M)=T(q)+KL(q:p). 



(3) 

(4) 
(5) 

(6) 



P (g\M) = If p{g,f,8\M)àfà8 



àfdd 



Ainsi F{q), appelée l'énergie libre de q par rapport à p, 
est une limite inférieure de p(g\A4) car KL(q : p) > 0. Par 
la suite, nous allons écrire l'expression de T{q) par 

T(q) = (lnp(g,f,8:M)) q + H(q) (7) 

où Tt(q) est l'entropie de q. [U H] 

2 Approche variât ionnelle 

Nous allons maintenant utiliser ces relation pour décrire 
le principe de l'approche variationnelle. L'idée de base est 
que l'utilisation directe de la loi a posteriori conjointe 
p{f \8\g; M) est souvent très coûteux pour, par exemple, 
être explorée par échantillonnage directe ou pour calculer 
les moyennes a posteriori f — J J f p(f, 8\g; A4) d8 df 



et 9 = j J 9p(f,9\g;M)dfd9. En effet, rare sont les 
cas où on puisse trouver des expressions analytiques pour 
ces intégrales. De même l'exploration de cette loi par des 
méthodes de Monté Carlo est aussi coûteuses. On cherche 
alors de l'approximer par une loi plus simple q(f , 9). Par 
simplicité, nous entendons par exemple une loi q qui soit 
séparable en / et en 9 : 



?(/,*) = <h (/)&(*) 



(8) 



Évidemment, cette approximation doit être fait de telle 
sorte qu'une mesure de distance entre q et p soit minimale. 
Si, d'une manière naturelle, on choisi KL(ç : p) comme 
cette mesure, on aura : 

(Qi,Q2) = argmin{KL(<2ig 2 : p)} = argmax{J r (ç 1 ç 2 )} 

(91,92) (91.92) 

(9) 

et sachant que KL (qiq2 '■ p) est convexe en qi à , q2 fixée 
et vise versa, on peut obtenir la solution d'une manière 
itérative : 

qi = argmin 9i {KL(gig 2 ■ p)} = argmax ?i {T(qiq 2 )} 
q 2 = argmin ?2 {KL(çiç 2 ■ p)} = argmax ?2 {T(qiq 2 )} 

(10) 

Utilisant la relation ([7]), il est facile de montrer que les 
solutions d'optimisation de de ces étapes sont 



92 (0) 



oc exp 
oc exp 



(lnp(g,f,d;M)} Î2ie) 
(lnp(g,f,0;M)) îi{f) 



(H) 



Une fois cet algorithme convergé vers qKf) et q^(9), on 
peut les utiliser d'une manière indépendante pour calculer, 
par exemple les moyennes f = J f 9Î(/) df et 

e* = je q* 2 {8)àe. 

Une deuxième étape de simplification est nécessaire pour 
être capable de calculer les espérances qui se trouvent dans 
ces exponentielles. Les calculs non paramétriques sont sou- 
vent trop coûteux. On choisit alors une forme paramé- 
trique pour ces lois de telle sorte qu'on puisse, à chaque 
itération, remettre à jours seulement les paramètres de ces 
lois, à condition cependant que ces formes ne changent pas 
au cours des itérations. La famille des lois exponentielles 
conjuguées ont cette propriété p] [H [H [6l [8] . Nous 
examinons ici, trois cas : 

2.1 Cas dégénérée 

Il s'agit de choisir pour qi(f) et q 2 (9) des formes dégé- 
nérées suivantes : 



Uf\f) = *(/-/) 

q 2 (9\9) =6(9-6) 



(12) 



Par conséquence, qu'au cours des itérations, nous aurons 
à remettre à jours / et 6 au cours des itérations, les pa- 
ramètres de la loi a posteriori jointe p(f, 9\g; A4). 
En remarquant alors que 



qi(f)<xp(9,f,9;M)<xp(f,0\g;M) 
q 2 (d)Kp(g,f,0;M)<xp(f,0\g;M) 



(13) 



que l'on utilise ensuite pour mettre à jour q 2 (9). On note 
alors que cet algorithme devient équivalent à ce qu'on peut 
apeller MAP Joint : 

/= argmaxj lp(g, f, 9; M) }= argmax^ lp(f,9\g;M) 

9= argmaxfl /, 9; M)j = argmaxj {p(f,9\g;M) 

(14) ' 

On remarque que l'on retrouve un algorithme du type 
MAP jointe. 

2.2 Cas particulier conduisant à l'algorithme 
EM 

Il s'agit de choisir, comme dans le cas précédent une 
forme dégénérée pour q 2 (9) = 5(9 — 9), ce qui donne 

ÏÏi(f) « p(g, f, 9; M) oc p(f, 9\g; M) oc p(f\9, g; M) 

. (15) 

ce qui signifie que qi(f) est une loi dans la même famille 
que la loi a posteriori p(f\g, 9; A4). Évidemment, si la 
forme de cette loi est simple, par exemple une gaussienne, 
(ce qui est le cas dans les situations que nous étudierons) 
les calculs seront simples. 

A chaque itération, on aurait alors à remettre à jours 9 
qui est ensuite utilisé pour trouver qi(f\9) = p(f\g, 9; A4), 
qui est utilisée pour calculer 



Q(9,9) = (lnp(g,f,9;M)l 



(16) 



Il est alors facile de voir que si p(f, 9\g; A4) est gaussienne 
à 9 fixé, on aura juste à calculer / = argmaxy? |p(/, 9\g; A4) 



9i(/|0) 

On remarque facilement l'équivalence avec l'algorithme 
EM. 

2.3 Cas particulier proposé pour les pro- 
blèmes inverses 

Il s'agit de choisir, pour qi(f) et q 2 (9) les mêmes fa- 
milles de lois que p(f\g, 9) et p(9\g, f), ce qui permet de 
profiter de la mise à jour facile de ces lois si des lois a 
priori correspondantes sont choisie dans les familles des 
lois conjuguées associé à la modélisation directe du pro- 
blème. 

Dans ce travail, dans un premier temps, nous allons 
considéré le cas des problèmes inverses linéaires : 

g = Hf + e (17) 

où H représente la forme discrétisé de la modélisation di- 
recte du problème et e représente l'ensemble des erreurs de 
mesure et de modélisation avec des hypothèses suivantes : 

p(g\H, f, 9 e ; M)=M(Hf, (l/6 e )I), 

p(f\ê f ;M)=Af(0,(l/9 f )(D t f D f )- 1 ), 

p(9 e ;M)=g(a eQ ,f3 eQ ), 
p(0 f ;M)=g(a fo ,P fo ) 

où 9 = (9 e = l/o-^,9f = l/a 2 j). On obtient alors faci- 
lement les expressions de p(g, f\9; A4), p(f\g, 9; A4) et 
P{9\g,f; A4) qui sont : 

p(g, f\H, e ;M)=Af(Hf, (l/0 e )I) Af(0 7 (l/e^D^Df)- 1 ) 
p(f\g,H,e f ;M)=Af(f,è), 
p(d e ;M)=g(a e ,0 e ), 
p(6 f ;M)=g(a f ,P f ) 
\ (19) 



(18) 



où les expressions de /, X, (a e , /3 e ) et (S/, /?/) seront don- 
nées en annexe. 

3 Application en restauration d'image 



Dans le cas de la restauration d'image où H a une 
structure particulière, et où l'opération Hf représente une 
convolution de l'image / avec la réponse impulsionnellc h, 
la partie difficile et coûteuse de ces calculs est celle du cal- 
cul de / qui peut se faire à l'aide de la Transformée de 
Fourier rapide. 

De même, l'approche peut très facilement être éten- 
due pour le cas de la restauration aveugle ou myope où 
on cherche à la fois d'estimer la réponse pulsionnelle h, 
l'image / et les hyper-paramètres 9. Pour établir l'expres- 
sions des différentes lois dans ce cas, nous notons que le 
problème directe, suivant que l'on s'intéresse à / (déconvo- 
lution) ouà/i (identification de la réponse impulsionnelle), 
peux s'écrire 







[** < F* < B > F > * -+ 








— 


S„*' < F >*< B > g 




(25) 


q(0 ei ) 


= 


Q(a el ,(3 ei ) avec 








= 


a e0 + M/2, 








= 


/3 e0 + l/2 < ee* > Ul 




(26) 




= 


Ç(af,f3f) avec 






a f 




a f0 + N/2, 






fit 




/3 /0 + l/2Tr{g t Q<// t 


>}. 


(27) 


q(a wj ) 




G{aj,bj) avec 






Ot w j 




a w0 + 1/2, 










(3 w0 + 1/2 < wl > . 




(28) 


où A = 


: \ 


',{a w j,j = !,•■•, N } et B = 


: diag{/3 ei ,i 


= !,■■■ 



On a ainsi l'expression des différentes composante de 
la loi séparable approchante. On peut en déduire facile- 
ment les moyennes des ces lois, car ces lois sont, soit des 
gaussiennes, soit des lois gamma. 



g(r)=h(r) * f(r) + e(r) = f(r) * h(r) + e(r) 
g=H f + e = F h + e 



(20) 



Pour permettre d'obtenir une solution bayésienne pour 
l'étape de l'identification, nous devons aussi modéliser h. 
Une solution est de supposer h = <&w où la matrice <I? 
est une matrice telle que &w représente la convolution 
(j)(r) *w(r). Ainsi les colonnes de <ï> représentent une base 
et les éléments du vecteur w représentent les coéfficients 
de la décomposition de h sur cette base. On a ainsi 

g{r)={(j) * w) * f(r) + e(r) = /*(</>* w)(r) + e(r) 
g=* Wf + e = F$w + e 

^ (21) 

Le problème de la déconvolution aveugle se ramène à l'es- 
timation de / et w avec des lois 



< w >=n w , 

< f >=/*/, 

< a f >= a f /Pf, 



< wlj >= 

< ff >■■ 



2 

wii 



[/*■ 

'2 



< (3 f >= a f /p 



(29) 



< ee* >= gg l - 2g[< F > * < w >]* 



(30) 

Pour le calcul des termes < W l W > et < F* F > qui in- 
terviennent dans les expressions de S/, S„ et [Fww t F t ] 
on peut utiliser le fait que F et W sont des matrices 
block-Toeplitz avec des blocs Toeplitz (TBT), on peut 
les approximer par des matrices block-circulantes avec des 
blocs circulantes (CBC) et les inverser en utilisant la TFD. 
Notons aussi que fij et n w peuvent être obtenu par opti- 
misation de 

J(jJi f )=\g - GW^fY <B>[g- <i>Wf] + (l/0/)||Q/|| 2 
J(nJ=[g - *FnJ <B>[g- <!>Fw] + \\Aw\\ 2 

(31) 



p(g\w, /, E £ ) = Af(*Wf, £ £ ) = Af(F<&w, £ £ ), 
avec £ £ = diag j j-, i = 1, • • • , M j et p(0 e i) = G(a e a, /3 e o) Les détails de ces calculs seront omis ici. 
p{f\0f)=M(O,{6 f D t f D f )- 1 ) a vecp(6 f )=g(a fo ,0fo), 
p(w\a) = 11^(0, avecp(a) = U j 0(a ,b )yj 

(22) 

Avec ces lois a priori , il est alors facile de trouver l'ex- 
pression de la loi conjointe p(f, w, e ,9f, a; g) et la loi a 
posteriori p(f,w,0 e ,6f,a\g). Cependant l'expression de 
cette loi 



Restauration avec Modélisation 
Gauss-Markov-Potts 



p(f,w,0 e ,6 f ,a\g) oc p(g\w, f,~E e )p(f\6 f )p(w\a) 
p(6 e )p(0 f )p(a) 

(23) 

n'est pas séparable en ses composantes. L'approche varia- 
tionnclle consiste donc à l'approximer par une loi sépa- 
rable 

p(f, w, 6 ei 6 f ,a\g) ~ q(f)q(w) Oj q( e ei)q{0 f ) Ylj q(otj) et 
avec les choix des lois a priori conjuguées en appliquant 
la procédure décrite plus haut, on obtient 



?(/) 
S/ 

/*/ 

q(w) 



= M {pif, S/) avec 

= [** < W* < B > W > *+ < f > Q*Q] \ 

= S/** <W >*< B > g, (24) 

= N(Hwi avec 



Le cas d'une modélisation gaussienne reste assez restric- 
tif pour la modélisation des images. Des modélisation par 
des champs de Markov composites (intensités-contours ou 
intensités-régions) sont mieux adaptées. Dans ce travail, 
nous examinons ce dernier. L'idée de base est de classer 
les pixels de l'images / = {f(r),r G 1Z} en K classes 
étiquettées par une variable discrète z(r) e {!,-■■ ,K}. 
L'image z(r) = {/(r),r G 1Z} représente ainsi la seg- 
mentation de l'image /(r). Chaque paquets des pixels 
fk = {/( r )i T * 6 T^-k } représente un ensemble fini des ré- 
gions compacts et disjointes : UilZki = T^k et UkTZk = T^- 
On suppose aussi que f k et fi^k^l sont indépendants. 

A chaque région est associée un contour. Si on repré- 
sente les contours de l'images par une variable binaire q(r), 
on a q(r) = à l'intérieure d'une région et q(r) = 1 aux 
frontières de ces régions. On note aussi que q(r) à partir 
de z(r) s'obtient d'une manière déterministe (voir Fig 1). 




FiG. 1 - Modèle de mélange et champs de Markov caché : 
image des intensités ou niveau de gris f(r), image z(r) 
de segmentation ou classification, image binaire q(r) des 
contours. 

Avec cette introduction, nous pouvons définir 

p(f(r)\z(r) = k, m k , v k ) = Af(m k ,v k ) (32) 

ce qui suggère un modèle de mélange de gaussienne pour 
les pixels de l'image 

P(.f( r )) = ^2 a kN(m k , v k ) avec a k = P(z(r) = k) (33) 
k 

Une première modélisation simple est donc, supposées que 
les pixels de l'images sont a priori indépendants, ce qui 
suggère 

p(z) =l[p(z(r)) (34) 
r 

Nous apcllons ce modèle, Mélange de Gaussiennes Indé- 
pendantes (MGI). 

Maintenant, pour prendre en compte la cohérence spa- 
tiale de ces pixels, nous devons introduire, d'une manière 
ou autre, une dépendance spatiale entre ces pixels. La mo- 
délisation markovienne est justement l'outil approprié. 

Cette dépendance spatiale peut être fait de trois ma- 
nières. Soit utiliser un modèle markovien pour z{r) et 
un modèle indépendant pour f(r)\z(r), soit un modèle 
markovien pour f(r)\z(r) et un modèle indépendant pour 
z(r), soit un modèle markovien pour f(r)\z(r) et un mo- 
dèle markovien aussi z(r). Nous avons examiné ces cas 
avec des modèles de Gauss-Markov pour f(r)\z(r) et le 
modèle de Potts pour z(r). Ce dernier peut s'écrire sous 
deux formes : 



p(z(r)\z(r') , r' £ V(r)) oc exp 



7 S(z(r)-z(r')) 
r'ev(r) 

(35) 



p(z) oc exp 7 ô(z(r) - z(r')) (36) 

renr<ev(r) 

Ces différents cas peuvent alors se résumer par : 
Modèle Gauss-Potts : 

p{f(r)\z(r) = k)= Af(m k ,v k ),Vr E K , , 
P(f\*)=Ilrenrt(m(r)Mr)) [ ' 

avec m(r) — m k ,Mr G lZ k et v(r) — v k ,Vr e lZ k , et p(z) 
Potts. 

Modèle de mélange indépendante de Gauss-Markov : 

p(f(r)\f(r'),r' e V(r),q(r,r')) = Af(m(r),v(r)),Vr e K 
p(f\z) oc Yl k Af(m k lk^k) 



avec lfe = l,Vr G 1Z k et une matrice de covariance. 
Modèle de Gauss-Markov-Potts : 

p(f(r)\f(r'),r> e V(r),ç(r,r')) = M(m(r),v(r)),W e K 
p(z(r)\z(r'),r' e V(r)) oc exp 7Er'eV(r) H z ( r ) - z ( r ')) 

(39) 

Quelque soit le modèle choisi parmi ces différents mo- 
dèles, l'objectif est d'estimer /, z et 6. Si on écrit la loi a 
posteriori jointe : 

p(g\f,o) P (f\z,e) P (z) 



p(f,z,6\g) = 



p(g\o) 



(40) 



P( z ) = UrP( z ( r ) = k ) = IL" 



T,ren S{z{r ^ mk) 



(38) 



et on cherche à l'approximer par une loi séparable g(/, z, 0\g) 
qi(f) q 2 {z) qs(e). 

Cependant, ici, nous choisissons d'approximer seulement 
p(f,z,0) par qi(f\z) qz{z) qs{8). Les détails de ces cal- 
culs et une comparaison de ces différents algorithmes sont 
en cours d'expérimentation et de rédaction et seront pu- 
bliés dans un délai très proche. 

5 Conclusion 

L'approche variationelle de l'approximation d'une loi 
par des lois séparables est appliquée au cas de l'estima- 
tion non supervisée des inconnues et des hyper-paramètres 
dans des problèmes inverses de restauration d'image (dé- 
convolution simple ou aveugle) avec des modélisations a 
priori gaussiennes, mélange de gaussiennes ou mélange de 
gaussiennes avec champ de labels markovien (champs de 
Markov caché). 
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